Bayesian Analysis of Solar Oscillations 
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ABSTRACT 

A Bayesian probability based approach is applied to the problem of detecting and parameter- 
izing oscillations in the upper solar atmosphere for the first time. Due to its statistical origin, this 
method provides a mechanism for determining the number of oscillations present, gives precise 
estimates of the oscillation parameters with a self-consistent statistical error analysis, and allows 
the oscillatory model signals to be reconstructed within these errors. 

A highly desirable feature of the Bayesian approach is the ability to resolve oscillations with 
extremely small frequency separations. The code is apphed to SOHO/CDS (Solar and Heho- 
spheric Observatory/Coronal Diagnostic Spectrometer) O V 629A observations and resolves four 
distinct P4, P5, Pg and P7 p- modes within the same sunspot transition region. This suggests 
that a spectrum of photospheric p-modes is able to propagate into the upper atmosphere of the 
Sun and Sun- like stars, and places precise observational constraints on models of umbral eigen 
modes. 



Subject headings: Methods: statistical 
oscillations — Waves 

Introduction 



Sun: helioseismology — Sun: oscillations — sunspots — Stars: 



The use of Bayesian methods has great poten- 
tial within astrophysics and has been ap plied in 
areas from binary stars to cosmology (see iLoredo 
[1990 , for an introduction) . Bayesian methods have 
recently been employed within some areas of solar 
physics, such as th e analysis of radi ochemical so- 
lar neutrino data (iSturrock fc Wheatland 2008) 



inver sion of Stokes profiles (jAsensio Ramos et al. 



[2007^ and an approach to solar flare prediction 
(iWheatlan d 2004 ). The work of IJavn es (1987), 
Bretthorsti ( 19881 ) and others has shown that the 
application of Bayesian statistical techniques to 
spectral analysis has many applications in physics, 
but it has not yet been exploited in solar physics 
research. Current analysis techniques applied to 
the problem of wave detection and parameteriza- 
tion in solar physics are not optimal to the prob- 
lem at hand. Particularly concerning the estima- 
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tion of oscillation parameters and their uncertain- 
ties, it is not clear how to interpret least squares 
fitting or the Fourier and wavelet transform with- 
out understanding a relation to probability theory. 
It is possible to extract much more information 
contained within the data by applying Bayesian 
statistical methods, compared to the traditional 
least squares, Fourier or wavelet analysis currently 
employed. The Bayesian method allows extremely 
precise estimates of the oscillation parameters to 
be made, with a consistent statistical analysis of 
their uncertainties. For example, the probability 
based approach allows the obtainable frequency 
resolution to be estimated, which is much higher 
than can be interpreted from a Fourier transform. 

We ap ply the methods described by Ijavned 
lil987l ) and lBretthors^ (|l988l ) to the problem of fre- 
quency estimation within solar data. A Bayesian 
numerical code is applied to artificial time se- 
ries data, typical of oscillations within the solar 
corona, to demonstrate the high precision param- 
eter estimation that can be achieved. It is shown 
that frequencies spaced closer than neighboring 



Fourier frequencies can be successfully resolved by 
a Bayesian model. This makes the Bayesian ap- 
proach ideal for determining the number of fre- 
quencies present in a time series. Section [S] applies 
the method to transition region data, demonstrat- 
ing that it is possible to detect and resolve the 
presence of multiple frequencies in a time series 
where a Fourier analysis is unable to do so and its 
application is invalid. 

2. Bayes' Theorem 

The posterior probability of a hypothesis H , 
given the data D and all other prior information 
/ is stated by Bayes' theorem: 



P{H\D,I) 



P{H\I)P{D\H,I) 



(1) 



Bayes' theorem derives from commutative logic 
and the prod uct rule of probability theory (see 
lGregoryll2005l ). Where P{H\T) is the prior proba- 
bility of H given /, or the prior; P{D\I) is the 
probability of the data given /, and is usually 
taken as a normalizing constant; P{D\H,I) is 
the direct probability of obtaining the data given 
the hypothesis and prior information. The direct 
probability is termed the sampling distribution, 
when the hypothesis is held constant and differ- 
ent sets of data are measured. This sampling dis- 
tribution has become the traditional approach to 
estimating the probability of oscillations within 
astrophysics, particularly within the field of so- 
lar physics. However, unlike a laboratory exper- 
imenter, or statistician, typically, we can obtain 
only one measurement of the process under obser- 
vation. To proceed, the current archetypal method 
is to assume that the data is one of a large num- 
ber of possible measurements from a given sam- 
pling space. This sampling space is estimated by 
the application of Monte-Carlo or Fisher-type ran- 
domization techniques to gene rate a large number 
of artificial 'datasets' (e.g . see.lO'Shea et al.ll2001 : 



Linnell Nemec &: Nemedll985l) . Assuming a par- 
ticular hypothesis, the probability of observing the 
data within this sampling space of artificial data 
is then used to estimate the level of confidence 
in the hypothesis. In the problem of oscillation 
detection, the level of confidence that there is no 
oscillating signal present within the data is usually 
estimated, the null hypothesis. 



Since we generally have only one measurement 
of the data, rather than generating a distribution 
of artificial 'observations', it appears more logical 
to test the probability of obtaining the measured 
data against different hypotheses, incorporating 
the prior information we have available. This is 
the basis of the Bayesian method. The direct 
probability is then termed the likelihood function 
when the data are considered constant and tested 
against different hypotheses. 

3. Application of Bayes theorem to oscil- 
lation detection 



Th is sect ion summa r izes t he results of iJavnes 
(|l987l ) and iBretthorstI (|l988l ) which are applied 
in the code to calculate the Bayesian probability 
density function and the results in the following 
sections. 

3.1. The likelihood function 

When applied to the question of oscillation de- 
tection, we wish to compute the probability of a 
particular time series model, given the data and 
all other prior information. To calculate the likeli- 
hood function, the probability of the noise must be 
calculated. If the true model was known, then the 
difference between the data and the model func- 
tion would be equal to the noise distribution. As- 
suming Gaussian distributed noise, the probability 
of obtaining a particular series of noise values e^ 
is given by: 



N 



P(ei...eAf|cr,/) oc J^ 



1 



V27 



-exp 



2^2 



(2) 
where N is the number of elements in the series, 
and a^ is the noise variance. The likelihood func- 
tion is then given by: 

1 ^ 
L{{B}, M, a) = a-^exp{-^Y.^d, fiU)?}, 

2 — 1 

(3) 
where di are the measured values of the data. We 
suppose that the measured data is a combination 
of the model function and the noise i.e. 

di = f{ti) +ei. 

Note that, the data sampling is not required to be 
evenly spaced, unlike the Fourier transform. 



In the most general case, the model as a func- 
tion of time can be expressed as: 



The orthonormal model functions are given by: 



/(i) = ^i?,G,(i,M), 



(4) 



i=i 



where Bj are the amplitudes, m is the total num- 
ber of component model functions Gj , which are 
functions of any number of other parameters {p}, 
such as frequency, decay rates. . . etc. Substituting 
for f{t), the summation in the likelihood (Eqn[3]) 
becomes: 






3=1 i=l 



J=l k=l 



(5) 

where the cross term for the general model func- 
tion /(t), in Eqn[31 can be expressed as a matrix of 
the component model function products summed 
over time gj^. i.e. 



N 



9jk 



Y,Gj{U)Gk{U), {l<j,k<m). (6) 



i=l 



The general model function f{t) in Eqn.[3]may 
be composed of any number of component model 
functions. This is more easily represented in ma- 
trix form by a square matrix with indices j ,k repre- 
senting the standard row-major matrix notation. 
Equation [5] is then greatly simplified if the matrix 
gjk is diagonal. 

3.2. Calculating orthonormal functions 

In the simplest case of an oscillating model 
function containing a single frequency /(i) = 
i3i cos(aji) -I- B2Sm{u!t), the matrix is essentially 
diagonal due to orthogonality. 

In a more complex model containing multiple 
oscillations, the matrix will not generally be diag- 
onal. To diagonalize the matrix, the component 
model functions in Eqn. [6] must be transformed 
to a set of orthogonal functions. The matrix gjk 
is always a symmetric m x m square matrix; any 
matrix of this form has m linearly independent 
orthonormal eigen vectors and is orthogonally di- 
agonalizable. 



-. 771 



(7) 



fc=i 



where Cjk is the fcth component of the jth normal- 
ized eigen vector of gjk with a corresponding eigen 
value Xj. The functions H{t) then satisfy the or- 
thonormality condition ^^^i Hj(ti)Hk(ti) = Sj^ 
where Sjk is the identity matrix. The general 
model equation (Eqn. [?]) can then be expressed 
using these orthonormal component model func- 
tions. The matrix of these functions is then diag- 
onal and Eqn. [5] is greatly simplified. 

3.3. Marginalized probability 

The marginalization process allows us to cal- 
culate the probability independently of the pa- 
rameters in which we may have no interest, such 
as the component model function amplitudes, 
noise. . . etc. Marginalization allows one to re- 
move parameters from further explicit consider- 
ation in the posterior distribution, by assigning 
prior probabilities and integrating the posterior 
probability distribution over the variable to be re- 
moved. The resulting marginal distribution has 
no explicit mention of the removed variable, but 
rather expresses the probability as a function of 
the remaining variables. IBretthorsti (jl988h de- 
rives the probability density as a function of the 
frequency parameters as follows. 

Expressing the summation from the likelihood 
function (Eqn. [5|) using the orthonormal model 
functions allows the likelihood function to be writ- 
ten in independent terms for each of the compo- 
nent model function amplitudes. The likelihood 
function is marginalized to be independent of the 
model amplitudes, by assigning a uniform prior 
and integrating over each of the amplitudes; this 
assumes that we have no prior information to con- 
strain the amplitudes of the component model 
functions. Assuming that we have no prior infor- 
mation to constrain the noise, it can be marginal- 
ized in a similar manner to the amplitudes by as- 
signing a Jeffreys prior and integrating over all 
positive values. These parameters are marginal- 
ized using uninformative priors, where they are 
not constrained to any particular values. This 
gives an upper limit to the uncertainty of the pa- 
rameter estimates. Should we have any prior in- 



formation to constrain the parameter prior proba- 
bilities, then greater precision estimates would be 
achieved. 

This process has a great advantage compared 
to least squared fitting, in that the probability is 
evaluated only as a function of the parameters of 
interest. Thus reducing the dimensionality of the 
computed parameter space, whereas all parame- 
ters must be considered simultaneously using a 
least squares approach. Even after marginaliza- 
tion of the posterior probability distribution, the 
Bayesian method still allows good estimates of the 
marginalized parameters to be recovered without 
intensive computation, as described in Sect.[4l 

3.4. The probability density function 

The resulting posterior probability density that 
a general oscillatory model is present within the 
data is given by: 



P{{lu}\D,I)^ 






(8) 



This probability density has been derived as a 
function of the angular frequency parameters only 
{w}, assuming data with an unknown noise vari- 
ance; where m is the number of component model 
functions, A'^ is the number of measurements in the 
data time series and d? is the mean square value 
of the data. 

It is the h"^ function which carries the frequency 
dependence of the probability density. Given by: 



h^ 



^ m 



where, 



N 



^j ^^diHjiU), 



(1 < J < m). 



(9) 



(10) 



The hj values are the projections of the data 
onto the orthonormal model functions defined by 
Eqn. [21 and h'^ is the mean square value of these 
projections as a function of {a;}. The maximum of 
this function gives the most probable frequency w, 
supported by the data, for each of the component 
functions assumed by the model. The correspond- 
ing maximum in the probability density function 



(Eqn. [S]) is sharply peaked at these frequency val- 
ues uj, since the form of the function is similar to 
an exponential. This allows very precise frequency 
estimates to be made, at a resolution much higher 
than can be estimated from the Fourier transform, 
as described in Sect. 14.2] and [5l 

In the simplest case where we assume a gen- 
eral model function containing a single stationary 
harmonic frequency given by f[t) = Bi cos(cji) -I- 
B2S\n{ujt)^ where sine and cosine are the compo- 
nent model functions Gj given in Eqn. [4l then 
the eigen values and eigen vectors of the matrix 
gjk described in Eqn. [7] are equal to \j = , 



and ejk = ±-^ respectively. The h'^ function is 
the exact general solution; if we approximate by 
neglecting the negligible non-diagonal elements of 



the matrix then /i^ 



N 



E 



N 



.j^i^^j-e*"* 1^ which 
is the Schuster periodogram. It is an important, 
but subtle, point that probability theory shows 
there is a direct relation between the Schuster pe- 
riodogram and the probability that there is a sin- 
gle harmo nic frequency w ithin the data. A s de- 
scribed bv I Javned (|l987t ): iBretthorstI (|l988l ). the 
maximum of the periodogram gives the most prob- 
able frequency assuming that: there is a single sta- 
tionary harmonic frequency present, the value of 
N is large, there is no constant component or low 
frequencie s, and the da t a has a white noise dis- 
tribution. [Ireland et al.l ( 20081 ) take advantage of 
this fact, by applying an algorithm for automated 
oscillation detection within the solar corona. 

4. Parameter Estimation 

Although, in the single frequency case, the pe- 
riodogram, Fourier transform and h'^ are similar; 
h^ has been derived from probability theory and 
can be understood in a statistical sense. This un- 
derstanding allows estimates of the model param- 
eters, and their precision, to be derived. These es- 
timates cannot be made from least squares fitting, 
the Fourier transform, or periodogram, alone with- 
out understanding their origin in probability the- 
ory. Even though the posterior probability density 
(Eqn. [5]) is derived to be independent of param- 
eters such as the noise variance and model am- 
plitudes, good estimates of these parameters and 
their uncertainties can be recovered due to the 
sharpness of the probability density around the 
most probable frequencies, as described in this sec- 



tion. The parameter uncertainties are not given 
directly by least squares fitting, or the Fourier 
transform, which would require a sampling distri- 
bution approach. As described in Sect. 2, this 
is computationally intensive, and it is question- 
able whether this approach is appropriate given 
a single measurement of the data. Here we out- 
line estimates of the model parameters and their 
variance. 

4.1. The expected noise variance (cr^) 

The Bayesian analysis allows the expectation 
value of the noise variance within the data to be 
calculated as: 



(-^> 



1 



N -m-2 



N 
i=l 






(11) 



4.2. The frequency parameters {uj} 

The most probable frequencies, of the applied 
model, are evaluated numerically from the loca- 
tion of the maximum within the probability den- 
sity function described in Sect. 13.41 The accuracy 
of the frequency parameters can be estimated by 
expanding h^ (Eqn. [9]) in a Taylor series. This ac- 
curacy is dependent on the Hessian matrix of h"^ 
evaluated at the most probable model frequencies 



Ojk 



2 dujj dujk ' 



{l<j,k<r). (13) 



The estimated angular frequency resolution is 
given by the variance of the probability density 
function for ojk: 



The expected noise variance is a function of 
CO and is estimated with the hj functions evalu- 
ated at the most probable frequencies to given by 
the probability density function. The expectation 
value of the noise variance is essentially the differ- 
ence between: the total square value of the data, 
and the total square value of the data projected 
onto the orthonormal model functions defined by 
Eqn. \7\ It is implicit in the Bayesian model that 
everything within the data that is not fitted by the 
model is assumed to be noise. Thus the expected 
noise variance gives an indication to what degree 
the model represents the data. We may increase 
the complexity of the model, by the addition of 
more component model functions. However, once 
the real signal within the data has been accounted 
for, the addition of more component functions will 
have the effect of reducing (cr^) by fitting the noise. 
A method to determine the point at which the 
model best represents the true signal is described 
in Sect.O 

The percentage accuracy e of the expected noise 
variance is given by: 



^^2/{N-m-A), 



(12) 



where N is the number of elements in the series 
and m is the number of component model func- 
tions. The standard deviation accuracy estimate 
of the expected noise variance is then equal to 



(OE— ' a<k<r) 

'-^ V-j 
j=l 3 



(14) 



where Ujk is the kth component of the jth eigen 
vector of bjk with a corresponding eigen value Wj, 
r is the number of model frequencies, and the ex- 
pected noise variance (tr^) is evaluated at the most 
probable frequencies. These most probable fre- 
quencies are then equal to tD ± Guj^. ■ 

Equations [T3] and [T3] show that the obtain- 
able frequency resolution is related to how sharply 
the probability density is peaked around the most 
probable frequencies and the magnitude of the 
noise variance within the data. As described in 
Sect. 13. 4[ the form of the probability density func- 
tion is sharply peaked around these frequencies; it 
is this sharpness which allows very high precision 
frequency estimates to be made, as described in 
Sect.[5l and permits the results obtained in Sect. [6] 

4.3. The amplitude parameters {B) 

If each oscillating function within the model 
is expressed in the form f{t) — Bcoscos{ujt) + 
Bsin sin(wi), then two component model functions 
(sine and cosine) are used to describe each fre- 
quency component in Eqn. 21 Where Bsin and 
Bcos represent the amplitude parameters of the 
sine and cosine functions. For multiple frequency 
models, 

f{t) = Bi cos(ajii) + B2 cos{uj2t) + 



Ba sin(ci;ii) + B4 sm{uj2t) . . . , 

where the Bk parameters of the cosine functions 
are indexed consecutively for each frequency com- 
ponent, followed by the sine functions. 

The expectation value of each amplitude pa- 
rameter is given by: 



Table 1: Obtained frequency resolutions for the 
single frequency data with a S/N=l. 

Analysis method Sf (mHz) 

Bayesian PDF (ct/) 002 

FFT (HWHM) 0.2 

Wavelet (HWHM) 0.5 



{Bk 



E 






(1 < fc < m), 



(15) 



where hj are the projections of the data onto the 
orthonormal model functions in Eqn. [7l Cjk are 
the components of the normalized eigen vectors of 
the matrix 5jfe given in Eqn.[Bl with corresponding 
eigen values Xj. These amplitude parameters are 
dependent on uj, and are estimated using the or- 
thonormal model functions evaluated at the most 
probable frequency parameters uj described in the 
previous section. This is a very good approxima- 
tion, due to the sharpness of the peak in the prob- 
ability density function which is almost described 
by a delta-function. 

The variance of each amplitude parameter is 
then given by: 



^k = 



N 



N-2 



2N-5 



2N-5~-2m 



2N ~7 



(P 



2N-7- 



2m 



N 



Ef^'(16) 



i=i 



where A'' is the number of time series elements, cP 
is the mean value of the data squared and /i^ is 
the mean square value of the data projected on to 
the orthonormal model functions. 

4.4. The polar amplitude (A) and phase 
parameters ((/>) 

The expected amplitudes can be used to express 
the model function results in polar coordinates, 
where each oscillation within the model is of the 
form: 

/(i) = A cos{ujt + (/)). 

The polar amplitude and phase for each frequency 
component are then given by: 



(A) = ./i32^, + i32.^, (<^) - arctan(^ 






■)• 



The Itj errors a a and a^ are then given by the 
propagation of the B^ Icr errors given in Eqn. 1161 



5. Frequency Resolution 

In this section, we apply the Bayesian model 
to artificial test data, typical of the type of oscil- 
lations that are observed within the solar corona 
using current instrumentation. 

5.1. The single frequency case 

Here we compare the results obtained from the 
Bayesian model with those obtained by perform- 
ing a Fourier, and wavelet analysis. The analyzed 
time series is of the form: 

di = cos{2TTfU) + sin{2TTf ti) + e^, 

where /=3.3 mHz and e^ is Gaussian distributed 
noise. This is typical of the 5-minute period 
band oscillations observed in coronal loops, with 
100 samples of 30 s cadence and has a low sig- 
nal to noise ratio with a RMS S/N=l. The 
Bayesian model is applied to calculate the prob- 
ability density function for a single harmonic fre- 
quency model of the time series, using Eqn.[8l The 
most probable frequency, given by the peak within 
the probability density function, can be described 
by a Gaussian of width a^j given by Eqn. [TH con- 
verting from angular frequency to give ct/ in Hz. 
This width determines the theoretical frequency 
resolution obtainable with the Bayesian model, 
which is much less than the estimated frequency 
resolution obtained using the Fourier transform. 

Figure [1] shows the obtainable resolving power 
of the Bayesian model, normalized to that ob- 
tained by the Fourier and wavelet transforms. The 
solid line indicates the probability density func- 
tion for the single frequency Bayesian model, the 
dashed line shows the FFT and the dotted line 
is the global wavelet transform using the Morlet 
wavelet function. We can see that the Bayesian 
model gives a significant increase in the resolu- 
tion of the estimated frequency, with the probabil- 
ity density almost described by a delta-function. 
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Fig. 1. — A comparison of the frequency resolu- 
tion obtained for the single frequency data with 
a S/N=l, with the normalized Bayesian probabil- 
ity density function (solid line) , FFT (dashed) and 
global wavelet transform (dotted). We see that the 
Bayesian PDF has an order of magnitude increase 
in resolution compared to the FFT. 

Table [T] lists the obtained frequency resolution 
for each method, estimating the resolution of the 
Fourier and global wavelet transforms using their 
half width half maximum (HWHM). We see that, 
for a S/N=l, the Itr error on the estimated fre- 
quency from the Bayesian model gives an order 
of magnitude increase in resolution over the FFT. 
The global wavelet has an even lower resolution 
due to the smoothing effect of the transform on 
the wavelet scale. In fact, if we are interested in 
high precision frequency measurements of station- 
ary frequencies, or closely separated frequencies, 
then a wavelet analysis is one of the worst meth- 
ods that we can apply. 

5.2. Two closely separated frequencies 

We now compare the results obtained from a 
Bayesian model and a Fourier analysis of two 
closely separated frequencies within a simulated 
time series shown in Fig [2^. Again, we generate 
a time series typical of coronal loop oscillations of 
the form: 

di = Bi cos{2TTfiti) + B2 cos(27r/2ii) -t- 
B3 s'm{2TTfiti) + B4 sm{2Trf2ti) + e^, 

with 100 samples at 30 s cadence, a harmonic fre- 
quency of /i =3.3 mHz, an additional frequency of 



Table 2: Estimated frequencies and la errors, of 
the parameters expressed in polar coordinates, ob- 
tained from the Bayesian model applied to the two 
frequency time series, separated by one Fourier fre- 
quency step and a RMS S/N=l. 



Frequencies f ±crf (mHz) 


True value (mHz) 


2.98 ± 0.05 
3.37 ± 0.05 


3.00 
3.33 


Amplitude A± aA 


True value 


0.96 ± 0.19 
0.94 ± 0.19 


1.00 
1.00 


Phase (/) ± 0-0 (rad) 


True value 


5.64 ± 0.14 
5.52 ± 0.15 


5.50 
5.50 


(-^> 


True value 


0.88 ± 0.13 


0.89 



/2=3.0 mHz, Gaussian distributed noise e^ and a 
RMS S/N=l. These two frequencies are separated 
by only one frequency step in the Fourier trans- 
form, so in principle their frequencies are directly 
adjacent in the FFT. 

Figure [21d shows the FFT for the two frequency 
time series; Fig. [5}; shows the Gaussian representa- 
tion of the resolution obtained from the Bayesian 
probability density function, which has been nor- 
malized to the FFT peak for comparison. Note 
Fig [21: is not a power spectrum, but an illustra- 
tion of the frequency resolution obtained with the 
Bayesian model. As expected the FFT is un- 
able to resolve two such closely separated frequen- 
cies. A single broad peak is observed, with a large 
HWHM, suggesting the possibility that more than 
one frequency may be present. However, the result 
from the Bayesian model resolves the two frequen- 
cies independently and to a very high precision 
even with a S/N=l. Table [2] lists the resolved 
frequencies and their la errors, estimated from 
the probability density function of the two har- 
monic frequency Bayesian model. We see that the 
Bayesian model not only resolves the two frequen- 
cies but docs so to a very high precision with Ict/ 
errors of 0.05 mHz, even with a relatively short 
duration time series. Figure [2jl shows the signal 
reconstructed from the Bayesian model parame- 
ters, and the true signal within the simulated time 
series. We see that the Bayesian code provides a 
very good reconstruction of the signal even with a 
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low signal to noise ratio. 

6. Application to Solar Data 

The Bayesian model is now ap plied to real ob- 



servat ions of solar oscillations. iMarsli fc Walsh 
( 20061 ) observe the apparent propagation of slow- 
magnetoacoustic waves within a sunspot region. 
These waves are observed to propagate from the 
transition region into the coronal loop system 
emerging from the sunspot and are interpreted as 
the propagation of photospheric p-niodes waveg- 
uided along the magnetic field. 

The original analysis applied Fourier tech- 
niques to the time series; here we apply the 
Bayesian model to th e O V data described in 



Marsh fc Walsh ( 



Marsh fc Walsh (2006). The results presented in 



2006( 1 show the presence of two 
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frequencies, in the 3-min period band, observed 
in the transition region above the sunspot um- 
bra. The data consist of 100 samples observed 
at 26 s cadence obtained with the Coronal Diag- 
nosti c Spectrometer (CDS) (see, [Harrison et al.l 
19951 ). The Bayesian model is applied iteratively, 
with the addition of further model functions to 
increase the complexity of the applied model. 

6.1. Model selection 

In addition to the problem of fitting a model 
to the data, we must determine which is the most 
probable model of those under test. If we already 
have a knowledge of the noise variance within the 
data, then the calculated expectation value of the 
noise variance (Eqnlll|l may be used to determine 
when the model has accounted for the real signal. 
Component functions may be added to the model 
until the expectation value of the noise variance 
equals the known noise variance. At this point 
the most appropriate model has been determined, 
and the addition of further component functions 
to the model will simply be fitting the noise. 

The noise pro perties of t h e CD S instrument 
are described by iThompsonI ([20001) • The noise 
variance of the CDS detector is given by ct^ = 
2Ap?i -|- R'^n, where Nph is the number of detected 
photons, R is the readout noise (here we use a con- 
servative value of 1 photon-event pixel"^), and n 
is the number of pixels summed over. We may use 
this known value of the noise variance to determine 
when the model has reached sufficient complexity 



Fig. 2. — a) Two frequency time series with a 
S/N=l. b) FFT of the time series containing two 
frequencies separated by 1 Fourier frequency step, 
c) Gaussian representations of the frequency reso- 
lution obtained from the Bayesian probability den- 
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Fig. 3. — Change in the expectation value of the 
noise variance for an increasing number of compo- 
nent frequencies within the Bayesian model. The 
dotted line indicates the known noise variance 
within the CDS data. 



able using a Fourier analysis. Fig. |4]d shows the 
Gaussian representation of the four frequencies re- 
solved by the Bayesian model. In addition to the 
oscillation frequency, the Bayesian code also re- 
turns high precision estimates of the amplitude 
and phase parameters, listed in Table [3l The 
height of the peaks in Fig. [4] are normalized to the 
corresponding amplitude parameters; note that 
this figure is not a power spectrum, but a rep- 
resentation of the frequency resolution obtained 
from the probability density function. As demon- 
strated in Sect. 15. 2|, the Bayesian model is able 
to resolve closely spaced frequencies to a much 
higher resolution than is possible with the Fourier 
transform, even with short duration observations. 
Where the FFT can resolve only two frequencies, 
the Bayesian model is able to resolve four indepen- 
dent frequencies within the data, to a very high 
precision. 



to account for the real signal and the addition of 
further model functions will begin to fit the noise. 
Figure [3] shows the change in the expectation 
value of the noise variance for increasingly com- 
plex models with the addition of more component 
frequencies. The error bars show the standard de- 
viation error estimate of the expected noise vari- 
ance derived from Eqn. [121 The dotted line shows 
the noise variance within the data due to the pho- 
ton statistics and the noise properties of the CDS 
detector. As expected, with the addition of fur- 
ther component frequencies to the model, the ex- 
pected noise variance is reduced. The expected 
noise variance reaches the level of the CDS data 
with a model containing four harmonic frequen- 
cies. Therefore we can state that the data best 
supports a model containing four frequencies. We 
are able to derive parameter estimates of these 
functions and their associated errors, from the 
probability density function of the four frequency 
model. 

6.2. The Bayesian results 

Figure [SK shows the FFT of t he O V data pre- 
sented in iMarsh &: WalshI ( 20061 ) , with two main 
peaks resolved in the transform. The broad width 
of the peaks may suggest that the data does not 
simply consist of two monochromatic frequencies, 
but this is the limit to the information avail- 



7. Conclusions 

When approaching the problem of oscillation 
detection, it is possible to extract much more 
information from the data by the application 
of a Bayesian model rather than the traditional 
least squares fitting, Fourier, or wavelet analysis. 
Considering the problem of frequency estimation 
within a time series, the Bayesian method returns 
very precise estimates and employs a rigorous 
self-consistent error analysis, due to its statisti- 
cal derivation. It is not clear how to determine 
frequency error estimates from the Fourier trans- 
form or least squares, without understanding their 
relation to probability theory. It is often miscon- 
ceived that the Fourier frequency spacing is the 
limit to the resolution with which a frequency can 
be resolved within a time scries. This is not the 
case, as the resolution limit is principally deter- 
mined by the signal to noise ratio. As shown in 
Sect. 15.11 we are able to estimate a single fre- 
quency, with S/N=l, to a resolution an order of 
magnitude greater than the FFT. In Sect. 15.21 we 
are able to resolve two oscillations with frequency 
separations directly adjacent in the FFT. This is 
not surprising as the Bayesian method has sim- 
ilarities with least squares fitting of a particular 
function. In an analogous way, the limiting resolu- 
tion with which an oscillation can be fitted with a 
sinusoidal function is not equal to the time series 
cadence, nor is the resolution with which a spec- 
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Fig. 4.— a) The FFT of the O V d ata originally 
presented in iMarsh fc WalshI ( 20061 1 showing two 
frequencies in the 3-niin range, b) Gaussian rep- 
resentations of the frequency resolution obtained 
from the Bayesian probability density function for 
the four frequency model, where the peaks are nor- 
malized to the derived amplitudes given in Table[3] 
c) Signal reconstructed from the Bayesian model 
(solid), original data (dot-dash). 



tral line can be fitted with a Gaussian equal to 
the pixel spacing on the detector; this resolution 
limit is largely determined by the signal to noise 
ratio within the data. 

describe how the 
results 



Loumos fc Deeminej (11975 

erroneous 



Fourier transform gives erroneous results for 
closely spaced frequencies. Their results a r e also 
explained by probability theory; iJavned ( 19871 ) 
discovered that the periodogram or Fourier trans- 
form is only directly related to the probability of 
a single stationary harmonic frequency within the 
data. If multiple frequencies are well separated, 
the Fourier transform still gives good frequency 
estimates, as the problem separates out into in- 
dependent single frequency probability problems. 
If the frequencies are closely spaced, however, the 
non-diagonal elements in the matrix gjk (Eqn. [6|) 
become significant and the frequencies are not 
orthogonal. The approximation of using a sin- 
gle frequency probability model for the purpose 
of frequency estimation is no longer valid, nor is 
the use of the Fourier transform. In this case, 
the transformation to orthogonal functions used 
in the Bayesian analysis is necessary to determine 
accurate frequency estimates and their uncertain- 
ties. 

As mentioned bv lBretthorstI ( 1988 ). the Bayesian 
method is similar to least squares, in that a least 
squares approach minimizes the summation in 
Eqn. [31 whereas the Bayesian results maximize 
the likelihood function. As described in Sect 13. 3[ 
the Bayesian method allows only the parameters of 
interest to be considered, greatly reducing the di- 
mensionality of the parameter space compared to 
a least squares approach. In the previous section, 
a least squares approach would require a 12 dimen- 
sional parameter space, where this is reduced to 4 
dimensions with the Bayesian model. In principal, 
least squares will produce similar frequency esti- 
mates to the Bayesian model, since uninformative 
priors have been used, however, least squares does 
not directly determine their uncertainties. The 
Bayesian framework allows the oscillatory model 
to be understood in terms of probability theory, 
and achieves higher precision estimates due to 
the sharp maximum of the posterior probability 
density. 

We apply the Bayesian model to the O V data 
presented in Marsh fc WalshI (|2006[ l. and resolve 
four closely spaced independent frequencies within 
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Table 3 

BAYESIAN model parameters derived from THe IMARSH fc WALShI (J2006l ) DATA. 



Frequencies / ± crj (niHz) 



Amplitude A ± a a (Photon-Events) Pliasc ± tr^ (rad) 



5.81 ± 0.05 
6.29 ± 0.03 
7.06 ± 0.05 
7.56 ± 0.04 



102.4 ± 26.2 
117.3 ± 26.1 
143.6 ± 26.2 
133.6 ± 26.3 



0.10 ± 0.18 
1.35 ± 0.16 
2.63 ± 0.13 
5.31 ± 0.14 



Expected noise variance (cr ) 
^ CDS noise variance a 



15062.0 ± 2271.0 
15029.0 



^CDS noise variance calculated using [Thonipsoni J2000I) 



the 3-min ute period range . The o bservations pre- 
sented in iMarsh &: WalshI (|2006[ ) are interpreted 
as the conversion and propagation of photosphcric 
p-mode oscillations along the magnetic field into 
the corona. 5-minute period range oscillations 
are observed within the umbral photosphere of 
sunspots, and are shown to be connected to the 
global p-m ode oscillation distribution centered on 
5-minutes ( Penn fc Labontelll993l : iBalthasar et al. 



which propagates through the umbral atmosphere. 
The frequencies detected here, and their spacing 
are consistent with the model results of Zhukov 



1987t iBraun et all Il987^ . Oscillations in the 3 
minute period range have been observed in the 
chrom osphere above sunspot umbrae for many- 
years (IBeckers fc Tallantill969 : Beckers fc Schultz 



19721 : iGurman et al.lll982l : rLites et al.lll982h . The 
3-minute period range oscillations are thought 
to be due to amplitude ste epening of th e pho- 
tosphcric p-i node spectrum ( Bogdanl 120001 ). Re- 
cent work bv lCenteno et alj ( 20061) supports this, 
demonstrating that 3-minute range power in the 
chromosphere is due to linear wave propagation 
from the 5-minute range power in the photosphere. 
The 3-minute oscillations a re also observed in 
the unibral transition region (IThomas et al.l 1987t 



Fludral l200ll lO'Shea et al] 120021: iRendtel et al 



20031 : iBrvnildsen et al.ll2004l ). The results pre- 
sented here suggest that we are able to resolve 
these oscilla tions into four closely spaced p-mode 
frequencies. IZhukovl (J2002l ) calculate the spectrum 
of eigen modes within the vertical magnetic field 
of the sunspot umbra, finding that the 3-minute 
umbral oscillations are due to p-modes modified 
by the s trong magnetic field within the sunspot. 
Zhukovl (J2005I ) calculates the same spectrum of 
umbral oscillations using the method of resonant 
filtering and by solving the eigen value problem, 
also determining that the 3-minutc oscillations 
are part of the photosphcric p-mode spectrum 



(|2005l ) and may represent the detection of the P4, 
P5, Pe and P7 photospheric p-modes in the solar 
transition region. These results provide precise 
observational constraints for future modeling of 
umbral eigen modes. A more detailed discussion, 
on the characterization of these resolved modes, 
is required than can be addressed here. A follow- 
ing paper will investigate these modes and their 
constraints on a model atmosphere. 
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